Ejercicio 1¶

1. Reporte descriptivo del dataset¶

Dimensiones del dataset, número de variables continuas y categóricas, distribuciones y comentarios generales

In [1]:
import pandas as pd
import numpy as np

pd.set_option('display.max_rows', None, 'display.max_columns', None)

fev = pd.read_csv(r'FEV_data.csv')

fev.head()
Out[1]:
seqnbr subjid age fev height sex smoke
0 1 301 9 1.708 57.0 2 2
1 2 451 8 1.724 67.5 2 2
2 3 501 7 1.720 54.5 2 2
3 4 642 9 1.558 53.0 1 2
4 5 901 9 1.895 57.0 1 2
In [2]:
fev.info()
fev.shape
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 654 entries, 0 to 653
Data columns (total 7 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   seqnbr  654 non-null    int64  
 1   subjid  654 non-null    int64  
 2   age     654 non-null    int64  
 3   fev     654 non-null    float64
 4   height  654 non-null    float64
 5   sex     654 non-null    int64  
 6   smoke   654 non-null    int64  
dtypes: float64(2), int64(5)
memory usage: 35.9 KB
Out[2]:
(654, 7)

Después de cargar el dataset y ver la información del mismo, vemos que tenemos 654 observaciones y 7 columnas, sin embargo la primera (seqnbr) es unicamente el número de secuencia o index de cada observación, mientras que subjid es el ID del paciente (en la buena teoría debería ser un número random). Por otro lado las columnas de sex y smoke son categóricas (realmente binarias), por lo que vamos a convertirlas al tipo 'category', no sin antes asegurarnos de que efectivamente contienen unicamente dos valores diferentes cada una.

In [3]:
fev.nunique()
Out[3]:
seqnbr    654
subjid    654
age        17
fev       575
height     56
sex         2
smoke       2
dtype: int64
In [4]:
factors = list(fev.loc[:, fev.nunique()<= 10])
fev[factors] = fev[factors].astype('category')
In [5]:
fev.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 654 entries, 0 to 653
Data columns (total 7 columns):
 #   Column  Non-Null Count  Dtype   
---  ------  --------------  -----   
 0   seqnbr  654 non-null    int64   
 1   subjid  654 non-null    int64   
 2   age     654 non-null    int64   
 3   fev     654 non-null    float64 
 4   height  654 non-null    float64 
 5   sex     654 non-null    category
 6   smoke   654 non-null    category
dtypes: category(2), float64(2), int64(3)
memory usage: 27.2 KB

Con esto, tendríamos una variable objetivo (fev) y 4 son posibles variables predictoras, de las cuales tenemos 2 variables continuas (age y heigth) y 2 variables categóricas.

No tenemos missing values puesto que todas las columnas tienen 654 datos, sin embargo debemos comprobar que no haya outliers, además queremos visualizar la distribución de las variables. Para esto haremos una tabla con datos descriptivos y posteriormente una inspección gráfica.

In [6]:
fev.describe()
Out[6]:
seqnbr subjid age fev height
count 654.00000 654.000000 654.000000 654.000000 654.000000
mean 327.50000 37169.571865 9.931193 2.636780 61.143578
std 188.93782 23690.860350 2.953935 0.867059 5.703513
min 1.00000 201.000000 3.000000 0.791000 46.000000
25% 164.25000 15811.000000 8.000000 1.981000 57.000000
50% 327.50000 36071.000000 10.000000 2.547500 61.500000
75% 490.75000 53638.500000 12.000000 3.118500 65.500000
max 654.00000 90001.000000 19.000000 5.793000 74.000000
In [7]:
fev.describe(exclude = np.number)
Out[7]:
sex smoke
count 654 654
unique 2 2
top 1 2
freq 336 589
In [8]:
import seaborn as sns
import matplotlib.pyplot as plt
import plotly.express as px

def histogram_boxplot(data, xlabel = None, title = None, font_scale=2, figsize=(9,8), bins = None):
    # Crear ventana para los subgráficos
    f2, (ax_box2, ax_hist2) = plt.subplots(2, sharex=True, gridspec_kw={"height_ratios": (.15, .85)}, figsize=figsize)
    # Crear boxplot
    sns.boxplot(x=data, ax=ax_box2)
    # Crear histograma
    sns.histplot(x=data, ax=ax_hist2, bins=bins) if bins else sns.histplot(x=data, ax=ax_hist2)
    # Pintar una línea con la media
    ax_hist2.axvline(np.mean(data),color='g',linestyle='-')
    # Pintar una línea con la mediana
    ax_hist2.axvline(np.median(data),color='y',linestyle='--')
    # Asignar título y nombre de eje si tal
    if xlabel: ax_hist2.set(xlabel=xlabel)
    if title: ax_box2.set(title=title, xlabel="")
    # Mostrar gráfico
    plt.show()


def plot(col):
     if col.dtypes != 'category':
        print('Variable Continua')
        histogram_boxplot(col, xlabel = col.name, title = 'Distibución continua')
     else:
        print('Cat')
        fig = px.bar(col.value_counts())
        fig.show()
        


fev.drop(['seqnbr', 'subjid'], axis = 1).apply(plot)
Variable Continua
Variable Continua
Variable Continua
Cat
Cat
Out[8]:
age       None
fev       None
height    None
sex       None
smoke     None
dtype: object

No tenemos missing data, no vemos errores en los datos; parece haber algunos outliers en la variable fev, sin embargo al tratarse de la variable objetivo y además, al no estar tampoco tan alejados de los demás valores al punto de parecer irreal, seguiremos trabajando con ellos para no perder posible información valiosa.

En cuanto a la distribución, tenemos para age y height distribuciones normales, fev, no es tan simétrica, pero de igual forma, posiblemente esté muy cerca de la normalidad. Mientas que las variables categóricas, tenemos por un lado Sex, que está muy bien distribuida entre sus dos categorías, con cantidad de datos muy similar para cada una; y por el otro lado, smoke (que como sería de esperarse al tratarse de niños y adolescentes) tiene muy poca representación de pacientes fumandores, la mayoría está en la categoría 2 (asumiremos que es no fumadores).

Parece ser una base de datos bastante limpia y lista para modelar!

2. Decidir si se descarta de inicio alguna de las variables de cara al modelado¶

Como se mencionó antes, la variable seqnbr, es el número de secuencia, por lo que no trabajaremos con ella.

La variable subjid pareciera ser también algun número de identificación del paciente, en ese caso no debería tomarse en cuenta puesto que se trata de números aleatorios.

La variable smoke, está muy poco representada en una de sus dos categorías, por lo que podría aportarnos muy poca información e incluso información sesgada, no vamos a descartarla desde ya, sino que trabajaremos con ella para ver si de alguna forma nos aporta información valiosa, pero está en la mira para descartarla en caso de que no nos aporte mucho.

3. Modelo de Regresión lineal sin interacciones entre las variables¶

Vamos a hacer una inspección visual y general, para hacernos una idea de las posibles relaciones entre las variables continuas

In [9]:
g = sns.PairGrid(fev.iloc[:, 2:])
g.map_upper(sns.scatterplot)
g.map_diag(sns.histplot)
g.map_lower(sns.regplot)
Out[9]:
<seaborn.axisgrid.PairGrid at 0x1c818948250>

Parece haber una clara relación entre las variables de edad y altura con la variable objetivo y también parece haber una dependencia entre ellas; sin embargo aún no hemos visto nada con respecto a sex y smoke que son categóricas.

Utilicemos pandas profiling para tener más información:

In [10]:
import pandas_profiling
pandas_profiling.ProfileReport(fev.iloc[:, 2:])
Summarize dataset:   0%|          | 0/5 [00:00<?, ?it/s]
Generate report structure:   0%|          | 0/1 [00:00<?, ?it/s]
Render HTML:   0%|          | 0/1 [00:00<?, ?it/s]
Out[10]:

En este reporte vemos resumido mucho de lo que analizamos antes, sin embargo hay datos adicionales en la sección de correlaciones: por un lado nos confirma una correlación fuerte entre age y heigth con la variable objetivo (y entre ellas también!) y por otro lado, nos indica que hay cierta relación entre sex y smoke con fev, aunque no es tan fuerte. Para esta segundo punto (correlación categórica-numérica) el reporte utiliza el método v de cramer; por lo que sería interesante ver este análisis en detalle.

Primero vamos a separar la variable objetivo para poder manipular las demás variables sin riesgo de modificar o transformar la variable de respuesta

In [11]:
varObj = fev.fev
inputs = fev.drop(['seqnbr', 'fev', 'subjid'],axis=1)
inputs.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 654 entries, 0 to 653
Data columns (total 4 columns):
 #   Column  Non-Null Count  Dtype   
---  ------  --------------  -----   
 0   age     654 non-null    int64   
 1   height  654 non-null    float64 
 2   sex     654 non-null    category
 3   smoke   654 non-null    category
dtypes: category(2), float64(1), int64(1)
memory usage: 11.9 KB

Y ahora si, usaremos la función de VCramer para un análisis preliminar de la relación de las variables con la variable objetivo (aunque ya tenemos una ligera idea de esto con base en la inspección visual que hicimos, pero queremos hacer la comparación de la correlación para cada input)

In [12]:
import scipy.stats as stats

def cramers_v(var1, vObj):
    
    if not var1.dtypes == 'category':
        bins = min(5,var1.value_counts().count())
        var1 = pd.cut(var1, bins = 5)
    if not vObj.dtypes == 'category':
        bins = min(5,vObj.value_counts().count())
        vObj = pd.cut(vObj, bins = 5)
        
    data = pd.crosstab(var1, vObj).values
    vCramer = stats.contingency.association(data, method = 'cramer')
    return vCramer


# Ejemplo uso univariante
#cramers_v(vinosCompra['Etiqueta'],vinosCompra['Beneficio'])

# Aplicar la función al input completo contra la objetivo
tablaCramer = pd.DataFrame(inputs.apply(lambda x: cramers_v(x,varObj)),columns=['VCramer'])

# Obtener el gráfico de importancia de las variables frente a la objetivo continua según vcramer
import plotly.express as px
px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones respecto al FEV').update_yaxes(categoryorder="total ascending")

Ahora si, podemos ver ordenados los inputs según su correlación con la variable objetivo!

A lo que vinimos, vamos a generar el modelo de regresión lineal para la predicción del FEV; lo primero es separar los datos en dos grupos: training y test; posteriormente generamos el modelo con los datos de training

In [13]:
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(inputs, varObj, test_size=0.2, random_state=42)

# Comprobamos dimensiones
print('Training dataset shape:', X_train.shape, y_train.shape)
print('Testing dataset shape:', X_test.shape, y_test.shape)
Training dataset shape: (523, 4) (523,)
Testing dataset shape: (131, 4) (131,)
In [14]:
# Unimos los datos de training de la variable objetivo, con los datos de training de los inputs:

fev_training = X_train.join(y_train)
fev_training.head()
Out[14]:
age height sex smoke fev
321 12 71.0 1 2 4.550
456 12 68.0 1 2 4.411
340 10 62.0 1 2 1.937
29 9 60.0 2 2 2.100
570 10 65.0 1 2 3.090
In [15]:
# Generamos el modelo completo:

from statsmodels.formula.api import ols 

reg_lineal = ols('fev ~ age + height + sex + smoke', data=fev_training).fit()
reg_lineal.summary()
Out[15]:
OLS Regression Results
Dep. Variable: fev R-squared: 0.766
Model: OLS Adj. R-squared: 0.764
Method: Least Squares F-statistic: 423.4
Date: Sun, 05 Mar 2023 Prob (F-statistic): 1.09e-161
Time: 09:45:19 Log-Likelihood: -281.13
No. Observations: 523 AIC: 572.3
Df Residuals: 518 BIC: 593.6
Df Model: 4
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -4.3576 0.271 -16.067 0.000 -4.890 -3.825
sex[T.2] -0.1511 0.038 -4.028 0.000 -0.225 -0.077
smoke[T.2] 0.0652 0.066 0.983 0.326 -0.065 0.196
age 0.0632 0.011 5.678 0.000 0.041 0.085
height 0.1041 0.005 19.150 0.000 0.093 0.115
Omnibus: 11.617 Durbin-Watson: 2.092
Prob(Omnibus): 0.003 Jarque-Bera (JB): 17.038
Skew: 0.178 Prob(JB): 0.000200
Kurtosis: 3.810 Cond. No. 932.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Como habíamos mencionado antes, la variable smoke está muy poco representada en una de sus categorías (es normal, los adolescentes no deberían fumar), además vemos que el p-value de su coeficiente de correlación es mayor a 0.05 (es el único mayor a 0); por lo que, como habíamos indicado antes, vamos a excluirla del modelo y ver cual es el efecto en el R cuadrado:

In [16]:
reg_lineal1 = ols('fev ~ age + height + sex', data=fev_training).fit()
reg_lineal1.summary()
Out[16]:
OLS Regression Results
Dep. Variable: fev R-squared: 0.765
Model: OLS Adj. R-squared: 0.764
Method: Least Squares F-statistic: 564.2
Date: Sun, 05 Mar 2023 Prob (F-statistic): 6.86e-163
Time: 09:45:52 Log-Likelihood: -281.61
No. Observations: 523 AIC: 571.2
Df Residuals: 519 BIC: 588.3
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -4.2829 0.260 -16.450 0.000 -4.794 -3.771
sex[T.2] -0.1540 0.037 -4.117 0.000 -0.227 -0.081
age 0.0601 0.011 5.629 0.000 0.039 0.081
height 0.1043 0.005 19.223 0.000 0.094 0.115
Omnibus: 11.573 Durbin-Watson: 2.100
Prob(Omnibus): 0.003 Jarque-Bera (JB): 17.840
Skew: 0.155 Prob(JB): 0.000134
Kurtosis: 3.850 Cond. No. 892.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

El R cuadrado del modelo practicamente no cambió al excluir smoke, por lo que lo mantendremos así. Parece ser que tenemos un buen modelo con un R2 de 0.765; sin embargo solo por entenderlo un poco mejor, vamos a ver cual es el aporte de cada variable al R2:

In [17]:
import patsy

# Generamos las matrices de diseño según la fórmula de modelo completo
y, X = patsy.dmatrices('fev ~ age + height + sex + smoke', fev_training, return_type='dataframe')
In [18]:
from relativeImp import relativeImp

names=X.columns.tolist()[1:]

# Calculamos importancia relativa
df_results = relativeImp(X.join(y), outcomeName = 'fev', driverNames = names)

# Ordenamos valores 
df_results.sort_values(by='normRelaImpt', ascending=False)
Out[18]:
driver rawRelaImpt normRelaImpt
3 height 0.443669 57.937958
2 age 0.277787 36.275651
0 sex[T.2] 0.024235 3.164859
1 smoke[T.2] 0.020075 2.621533
In [19]:
px.bar(df_results,x='normRelaImpt',y='driver',title='Importancia relativa por aportación al R2').update_yaxes(categoryorder="total ascending").show()

Según esto, podríamos eliminar (además de smoke) también a sex del modelo, veamos como cambia:

In [20]:
reg_lineal2 = ols('fev ~ age + height', data=fev_training).fit()
reg_lineal2.summary()
Out[20]:
OLS Regression Results
Dep. Variable: fev R-squared: 0.758
Model: OLS Adj. R-squared: 0.757
Method: Least Squares F-statistic: 812.9
Date: Sun, 05 Mar 2023 Prob (F-statistic): 8.87e-161
Time: 09:46:25 Log-Likelihood: -290.02
No. Observations: 523 AIC: 586.0
Df Residuals: 520 BIC: 598.8
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -4.5960 0.253 -18.181 0.000 -5.093 -4.099
age 0.0523 0.011 4.900 0.000 0.031 0.073
height 0.1095 0.005 20.417 0.000 0.099 0.120
Omnibus: 18.099 Durbin-Watson: 2.110
Prob(Omnibus): 0.000 Jarque-Bera (JB): 31.130
Skew: 0.231 Prob(JB): 1.74e-07
Kurtosis: 4.103 Cond. No. 852.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

El R2 se disminuye ligeramente, pero nada considerable, por lo que tenemos dos posibles modelos lineales: uno que incluye edad, altura y sexo, y otro únicamente con edad y altura (que desde hace rato vimos que eran las que parecían tener mayor correlación).

Vamos a evaluar ambos modelos con los datos de test:

In [21]:
from sklearn.metrics import mean_squared_error, r2_score

# Predicciones para test modelo completo
fev_pred1 = reg_lineal1.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred1)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred1))
Mean squared error: 0.40
Coefficient of determination: 0.80
In [22]:
# Predicciones para test modelo completo
fev_pred2 = reg_lineal2.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred2)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred2))
Mean squared error: 0.41
Coefficient of determination: 0.79

Tenemos datos muy similares para ambos modelos; los dos son aceptables, pero una variable menos nos simplificaría mucho el modelo. Veamos las interacciones:

4. Ajusta el mejor modelo de regresión lineal con interacciones entre las variables que resulten relevantes.¶

Ya habíamos visto en nuestro análisis gráfico, que hay una aparente interacción entre age y height, que además como vimos antes, vienen siendo las variables más relevantes para la predicción del valor de FEV, por lo que vamos a tomar en cuenta esa interacción:

In [23]:
reg_lin_interac2 = ols('fev ~ age + height + age*height', data= fev_training).fit()
reg_lin_interac2.summary()
Out[23]:
OLS Regression Results
Dep. Variable: fev R-squared: 0.783
Model: OLS Adj. R-squared: 0.781
Method: Least Squares F-statistic: 622.7
Date: Sun, 05 Mar 2023 Prob (F-statistic): 1.72e-171
Time: 09:48:36 Log-Likelihood: -261.65
No. Observations: 523 AIC: 531.3
Df Residuals: 519 BIC: 548.3
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -0.4138 0.593 -0.698 0.485 -1.579 0.751
age -0.4306 0.063 -6.789 0.000 -0.555 -0.306
height 0.0410 0.010 4.013 0.000 0.021 0.061
age:height 0.0077 0.001 7.712 0.000 0.006 0.010
Omnibus: 17.104 Durbin-Watson: 2.071
Prob(Omnibus): 0.000 Jarque-Bera (JB): 37.858
Skew: -0.053 Prob(JB): 6.02e-09
Kurtosis: 4.314 Cond. No. 2.25e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.25e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
In [24]:
reg_lin_interac1 = ols('fev ~ age + height + sex + age*height', data= fev_training).fit()
reg_lin_interac1.summary()
Out[24]:
OLS Regression Results
Dep. Variable: fev R-squared: 0.785
Model: OLS Adj. R-squared: 0.784
Method: Least Squares F-statistic: 474.1
Date: Sun, 05 Mar 2023 Prob (F-statistic): 1.49e-171
Time: 09:48:40 Log-Likelihood: -258.17
No. Observations: 523 AIC: 526.3
Df Residuals: 518 BIC: 547.6
Df Model: 4
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -0.5429 0.592 -0.918 0.359 -1.705 0.619
sex[T.2] -0.0968 0.037 -2.635 0.009 -0.169 -0.025
age -0.3881 0.065 -5.961 0.000 -0.516 -0.260
height 0.0431 0.010 4.229 0.000 0.023 0.063
age:height 0.0071 0.001 6.971 0.000 0.005 0.009
Omnibus: 14.277 Durbin-Watson: 2.066
Prob(Omnibus): 0.001 Jarque-Bera (JB): 28.474
Skew: -0.054 Prob(JB): 6.56e-07
Kurtosis: 4.138 Cond. No. 2.26e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.26e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

De nuevo vemos que la presencia o no de sex en el modelo, no afecta mucho, de hecho afecta aún menos en el modelo que incluye la interacción de edad y altura (que por cierto, parece ser un modelo levemente mejor que el que no la incluye si nos basamos en el R2.

Vamos a ver la evaluación de estos modelos al comparar las predicciones con los datos de test y ya luego veremos la validación cruzada de ambos modelos:

In [25]:
# Predicciones para test modelo completo
fev_pred3 = reg_lin_interac1.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred3)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred3))
Mean squared error: 0.39
Coefficient of determination: 0.82
In [26]:
# Predicciones para test modelo completo
fev_pred4 = reg_lin_interac2.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred4)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred4))
Mean squared error: 0.39
Coefficient of determination: 0.81

Efectivamente vemos que el aporte de sex al modelo, es mínimo, podemos excluirlo para simplificar el modelo sin perder información importante!

Comparación de ambos modelos por validación cruzada repetida. ¿Cuál de ellos tiene mejor comportamiento en generalización?¶

In [27]:
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from sklearn.linear_model import LinearRegression

def cross_val_lin(formula, data, seed=12345):
    y, X = patsy.dmatrices(formula, data, return_type='dataframe')
    model = LinearRegression()
    # Establecemos esquema de validación fijando random_state (reproducibilidad)
    cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=seed)
  
    # Obtenemos los resultados de R2 para cada partición tr-tst
    scores = cross_val_score(model, X, y, cv=cv)
  
    # Sesgo y varianza
    print('Modelo: ' + formula)
    print('Coeficiente de determinación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
      
    #sns.violinplot(y=scores,palette='viridis')
      
    return(scores)
In [28]:
form1 = 'fev ~ age + height'
form2 = 'fev ~ age + height + age*height'



# Creamos lista de fórmulas   
list_form = [form1,form2]
list_form

# Aplicamos a toda la lista la función creada (devuelve un dataframe pero está transpuesto)
list_res = pd.DataFrame(map(lambda x: cross_val_lin(x,fev, seed=2022),list_form))
list_res
Modelo: fev ~ age + height
Coeficiente de determinación R2: 0.760 (0.030)
Modelo: fev ~ age + height + age*height
Coeficiente de determinación R2: 0.783 (0.033)
Out[28]:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
0 0.787667 0.716968 0.799950 0.729255 0.761944 0.791089 0.776350 0.688730 0.788862 0.765495 0.772354 0.786487 0.769036 0.718342 0.762030 0.784979 0.741495 0.785038 0.709679 0.772243 0.750218 0.789058 0.750301 0.784801 0.719166 0.711838 0.761917 0.725076 0.818966 0.767828 0.723102 0.712695 0.774206 0.809178 0.785401 0.823126 0.691026 0.811060 0.708003 0.747786 0.797421 0.736606 0.727761 0.768323 0.780752 0.777619 0.757366 0.727513 0.764900 0.772503 0.802393 0.762426 0.778155 0.698572 0.755787 0.803356 0.710792 0.782174 0.763558 0.727841 0.809943 0.768817 0.756976 0.748240 0.746915 0.73496 0.757857 0.791671 0.743574 0.768408 0.794476 0.765307 0.778449 0.736414 0.724770 0.760133 0.738737 0.783343 0.806688 0.720035 0.726006 0.721742 0.775253 0.781166 0.775853 0.731595 0.778746 0.72327 0.784665 0.798909 0.771248 0.750733 0.771886 0.727637 0.769702 0.775609 0.761466 0.763527 0.750511 0.760841
1 0.808889 0.743159 0.815954 0.749891 0.798461 0.811031 0.812616 0.721601 0.811489 0.762009 0.801723 0.817343 0.795215 0.740051 0.777243 0.813767 0.770766 0.816728 0.697378 0.800250 0.794001 0.801966 0.779349 0.807909 0.730690 0.716065 0.801887 0.733507 0.845703 0.796791 0.737427 0.746988 0.790016 0.819267 0.825149 0.840542 0.706053 0.840481 0.726434 0.780371 0.834036 0.756675 0.741305 0.791445 0.800239 0.813321 0.780701 0.724106 0.796349 0.794709 0.810594 0.769931 0.805100 0.734141 0.790558 0.836799 0.744945 0.803868 0.771287 0.752825 0.822430 0.779777 0.786041 0.788226 0.765887 0.75224 0.770672 0.815152 0.788136 0.789405 0.821415 0.789580 0.812760 0.758850 0.733953 0.779442 0.765022 0.805122 0.824413 0.753467 0.731746 0.748743 0.796179 0.817650 0.799745 0.743433 0.785447 0.77693 0.823026 0.790566 0.796992 0.751288 0.806779 0.759025 0.791933 0.797006 0.793718 0.756039 0.793153 0.785192
In [29]:
results = list_res.T.melt()
results.columns = ['Modelo','R2']
results.groupby(['Modelo']).describe().T
Out[29]:
Modelo 0 1
R2 count 100.000000 100.000000
mean 0.760046 0.782857
std 0.030248 0.032920
min 0.688730 0.697378
25% 0.736050 0.756516
50% 0.764229 0.790287
75% 0.781418 0.807061
max 0.823126 0.845703

Vemos un mejor coeficiente de determinación para el modelo con interacciones, aunque también una leve desmejora en el error estandar, pero esta desmejora es mínima, entonces nos vamos a dejar influenciar más por el coeficiente de determinación (que tampoco tiene diferencia abismal, pero si se ve mejoría). Veamos de una forma gráfica la diferencia en el R2:

In [30]:
sns.boxplot(x='Modelo',y='R2',data=results,palette='viridis')
Out[30]:
<AxesSubplot:xlabel='Modelo', ylabel='R2'>

Por lo que podemos ver, ambos modelos son relativamaente similares, aunque el modelo con form2 (con interacción) es el de mejor ajuste por lo que este sería nuestro modelo final!

form2 = 'fev ~ age + height + age*height'